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Abstract 

In this paper, we study the synthesis of Gegenbauer processes using the wavelet 
packets transform. In order to simulate a 1-factor Gegenbauer process, we introduce 
an original algorithm, inspired by the one proposed by Coifman and Wickerhauser 
[1], to adaptively search for the best-ortho-basis in the wavelet packet library where 
the covariance matrix of the transformed process is nearly diagonal. Our method 
clearly outperforms the one recently proposed by [2], is very fast, does not depend 
on the wavelet choice, and is not very sensitive to the length of the time series. 
From these first results we propose an algorithm to build bases to simulate /c-factor 
Gegenbauer processes. Given its practical simplicity, we feel the general practitioner 
will be attracted to our simulator. Finally we evaluate the approximation due to the 
fact that we consider the wavelet packet coefficients as uncorrelated. An empirical 
study is carried out which supports our results. 
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1 Introduction 



The simulation of long memory processes is an issue of a paramount impor- 
tance in many statistical problems. In the time domain, there exist different 
methods devoted to this task (see [3] for a non exhaustive review of them). 
Alternative efficient approaches, which operate in the frequency domain, were 
also proposed (see [4,5,3]). More recently, owing to their scale-invariance prop- 
erty, wavelets have since been widely adopted as a natural tool for analyzing 
and synthesizing 1/f long-memory processes. They were demonstrated to pro- 
vide almost Karhunen-Loeve expansion of such processes [6]. 

The simulation of fractional differenced Gaussian noise (fdGn) using discrete 
wavelet transform (DWT) has been studied by [7]. This kind of process is 
characterized by an unbounded power spectral density (PSD) at zero. The 
proposed method relies on the fact that the DWT approximately decorre- 
lates long memory processes (see e.g. [8,9,6,10,11]). The orthonormal wavelet 
decomposition "only" ensures approximate decorrelation. The quality of this 
approximation has been widely assessed in [12,13,14,9,6,10,15] for a variety of 
1/f long memory processes. 

The DWT is only adapted to processes whose PSD is unbounded at the origin. 
Gegenbauer processes (sometimes also called seasonal persistent processes) are 
also long memory processes and are characterized by an unbounded PSD. The 
main difference with the fdGn processes is that the singularities of the PSD 
of the Gegenbauer processes can be located at one or many frequencies in the 
Nyquist domain, not necessarily at the origin. Therefore, a natural tool to an- 
alyze such processes appears to be the wavelet packet transform (WPT) , which 
is a generalization of the wavelet transform. The wavelets packets adaptively 
divide the frequency axis into separate dyadic intervals of various sizes. They 
segment unconditionally, the frequency axis and are uniformly translated in 
time. Moreover, a discrete time series of size N is decomposed into more than 
2 Ar / 2 wavelet packet (WP) bases. Among these bases, one is a very good can- 
didate to whiten the series and hence almost diagonalizes the covariance of 
the seasonal process. 

Recently, Mallat, Zhang and Papanicolaou [16], and, following their work, 
Donoho, Mallat and von Sachs [17], studied the idea of estimating the covari- 
ance of locally stationary processes by approximating the covariance of the 
process by a covariance which is almost diagonal in a specially constructed 
basis (cosine packets for their locally stationary processes) using an adap- 
tation of Coifman-Wickerhauser (CW) best ortho-basis algorithm. To some 
extent (given that we are interested in synthesis and they were in estimation 
issues), our work here can be seen as the spectral dual of theirs, since we are 
interested in studying the covariance of seasonal processes in the WP domain. 
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To the best of our knowledge, the simulation of the Gegenbauer process using 
the Discrete WPT (DWPT) has been first studied in [2]. The DWPT creates a 
redundant collection of wavelet coefficients at each level of the transform orga- 
nized in a binary tree structure, equipped with a natural inheritance property. 
Different methods exist to determine the best candidate orthonormal basis. 
The author in [2] used a method which depends on both the location of the 
singularity and the wavelet used in the DWPT. To simulate realizations of 
a Gegenbauer process, once the basis is found, it then remains to apply the 
(inverse) DWPT using the same approximation as in [7]. 

This basis search method consists first in considering the square gain function 
of the wavelet filter associated with each WP coefficient that is sufficiently 
small at the Gegenbauer frequency. Then a pruning of this family is done to 
obtain the ortho-basis. The main advantage of this method is its simplicity. 
However several points are still questionable and must be clarified. First the 
notion "sufficiently small" implies the introduction of a threshold which seems 
to depend both on the wavelet used and the length of the simulated series. 
No indication is given how to choose this threshold which remains awkward 
to control. Furthermore, it is not clear why the basis should depend on the 
wavelet. Lastly, this method inherently leads to an over-partitioning of the 
spectra which depends on the wavelet and the threshold considered (see e.g. 
Fig. 2 and Fig.4 in [2]). Indeed, as the Gegenbauer process we consider here 
is stationary, it is known that the Karhunen-Loeve basis is the Fourier ba- 
sis. While over-partitioning, the approach of [2] inherently tries to approach 
the Fourier basis (more precisely it tends to select most of the atoms from the 
Shannon wavelet packets at the deepest level). Then, this makes wavelet pack- 
ets machinery only of limited interest here. Furthermore, many important sta- 
tistical tasks involving Gegenbauer processes would seriously suffer from such 
an over-partitioning, e.g. maximum likelihood estimation, resampling-based 
inference, to cite only a few examples. 

To alleviate these intricacies, our belief is that it should more beneficial to 
build, for each Gegenbauer process, an unique valid (almost whitening) ba- 
sis for all wavelets with a reduced number of packets. This basis should only 
depend on the Gegenbauer frequencies, but not on the long memory param- 
eters nor on the wavelet used. The rationale behind these claims can be sup- 
ported by different arguments. Indeed, wavelets are now widespread as almost- 
diagonalizing expansion for 1/f processes, no matter what the long memory 
parameter and the wavelet are. Although, the latter parameters clearly in- 
fluence the quality of the decorrelation as was widely proven [6]. Our goal is 
then to mimic this behavior by extending and generalizing the aforementioned 
properties to Gegenbauer processes within the WP framework, with the desir- 
able properties that (i) our basis tends to the dyadic wavelet basis (the 1-band 
WP) when the singularity frequency tends to 0, and (ii) the provided basis 
should have a limited number of packets. To get a gist of the latter property, 
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we can say that we are seeking a basis (and the corresponding tree) which 
attains the minimal diagonalization error penalized by the complexity of the 
tree in terms of the number of packets (i.e. number of leaves of the tree) in- 
volved in the dyadic partition of the spectral axis provided by the selected WP 
basis. See [17, Sec. 13] and [18] for a more detailed discussion of complexity 
penalized estimation and its relation to best-ortho-basis. 

In this paper, we propose an alternative efficient way to determine the appro- 
priate basis for the simulation of 1-factor Gegenbauer process, that we extend 
to the simulation of /c-factor Gegenbauer process. To find this basis, We pro- 
pose an algorithm which is an adaptation of the best-basis search algorithm of 
[1]. The main property of this algorithm is that it provides us with a (unique) 
basis that only depends on the location of the Gegenbauer frequency, unlike 
the construction method of [2] which provides bases depending both on the 
location of the singularity and the wavelet. To point out the role played by 
the wavelet used, we will study the decorrelation properties of the WP coeffi- 
cients of a Gegenbauer process when it is expressed in this basis. In particular, 
the influence of the wavelet regularity, the long memory parameter and the 
location of the singularity on the decorrelation decay speed will be established. 

The organization of this paper is as follows. After some preliminaries and nota- 
tions related to the WPT theory (Section 2.1) and to the Gegenbauer process 
(Section 2.2) are introduced, we will define the best-basis search algorithm 
and the cost function we propose (Sections 3.1-3.2). Theoretical support to 
this cost function is also supplied. We then develop an algorithm to build an 
appropriate basis to simulate 1-factor Gegenbauer process (Section 3.3). This 
method will then be extended to fc-factor processes (Section 3.4). Theoretical 
evaluation of the approximation quality due to the fact that we consider the 
WP coefficient as uncorrelated is studied in Section 4. A simulation study is 
finally conducted to illustrate and discuss our results (Section 5). 



2 Preliminaries 



2. 1 The wavelet packet transform 



Wavelet packets were introduced by Coifman, Meyer and Wickerhauser [19], by 
generalizing the link between multi-resolution approximations and wavelets. 
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Let the sequence of functions defined recursively as follows: 

oo 

V&i(*)= E KnW^t-Vn) (1) 

n=— oo 

oo 

E g{n)Wt-#n) (2) 

n=— oo 

for j 6 N and p — 0, . . . , 2 J ' — 1, where /z and g are the conjugate pair of 
quadrature mirror filters (QMF). At the first scale, the functions ipo and f/'i 
can be respectively identified with the father and the mother wavelets <fi and 
tpi with the classical properties (among others): 

j<p(t) = i,Jm = o (3) 

The collection of translated, dilated and normalized functions tpj' n = f 2~i/ 2 ip p {2~H— 
n) makes up what we call the (multi-scale) wavelet packets associated to the 
QMFs h and g. j G N is the scale index, p = 0, . . . , 2 3 ' — 1 can be identified 
with a frequency index and k is the position index. It has been proved (see 
e.g. [20]) that if {ipf n } n( zz is an orthonormal basis of a space Vj, then the 
family {if)j+i, ^f+i^lnez is also an orthonormal basis of Vj. 

The recursive splitting of vector spaces is represented in a binary tree. To each 
node (j,p), with j G N and p = 0, . . . , 2 J — 1, we associate a space with 
the orthonormal basis {ip^{t — 2^n)} ne z- As the splitting relations creates two 
orthogonal basis, it is obvious that = © V^ 1 . 

The WP representation is overcomplete. That is, there are many subsets of 
wavelet packets which constitute orthonormal bases for the original space Vo 
(typically more than 2 2 for a binary tree of depth J). While they form a 
large library, these bases can be easily organized in a binary tree and efficiently 
searched for extreme points of certain cost functions, see [1] for details. Such 
a search algorithm and associated cost function are at the heart of this paper. 

In the following we call the collection B = {^' n }^^) € r,nez the basis of L 2 (R), 
and the tree T for which the collection of nodes (j,p) are the leaves, the 
associated tree. 

Given a basis B and its associated tree T it is possible to decompose any 
function x of L 2 (M) in B. At each node (j,p) G T, the WP coefficients Wj{n) 
of x in the subspace at position n are given by the inner product: 

Wj{n) = J i)){t - 2 j n)x(t)dt. (4) 
For a discrete signal of iV equally-spaced samples, the DWPT is calculated 
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using a fast filter bank algorithm that requires 0(N log N) operations. The 
interested reader may refer to the books of Mallat [21] and Wickerhauser [20] 
for more details about the DWPT. 

2.2 Gegenbauer process 

The fc-factor Gegenbauer process is a l//-type process introduced in [22,23]. 
The PSD / of a such process (X t )t is given by for all |A| < 1/2 

2 k 

= y n ( 4 ( c ° s 2?rA - c ° s 2 ™tf) ~* ( 5 ) 

i=l 

where k is a finite integer and 0<rfj<l/2if0< | i/j | < 1/2 and < dj < 1/4 
if |z/j| = for z = 1, . . . , fc. The parameter di and z/j are respectively called the 
memory parameter and the Gegenbauer frequency. The fc-factor Gegenbauer 
process is a generalization of the fractionally differenced Gaussian white noise 
process (see [24] and [25]) in the sense that the PSD is unbounded at k different 
frequencies not necessary located in 0. 

The Gegenbauer process (X t )t is related to a white noise process (e t )t with 
mean and variance a 2 through the relationship: 

f[(I-2v i B + B 2 ) d *X t = £t , (6) 

i=l 

where BX t = X t -± and rji = cos27rz/j. 

The main characteristic of the Gegenbauer processes in the time domain is the 
slow decay of autocovariance function. In the case of a 1-factor Gegenbauer 
process, Gray et al. [22] and then Chung [26] proved the asymptotic behavior 
of the autocovariance function: 

p(h) ~ h 2d ~ x cos(2iruh) as h -»• oo. (7) 

The next section is devoted to the construction of the best basis diagonalizing 
the covariance of a iV-sample realization of a Gegenbauer process with the 
convention N = 2 J . 



3 Simulation of Gegenbauer processes 

This section is composed of two parts. The first one is devoted to the simu- 
lation procedure in the general case: no assumption is made concerning the 
basis, except that we have an appropriate basis. The second part concerns the 
construction of this appropriate basis. 
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3. 1 Simulation procedure 

Here we present the procedure to simulate a Gegenbauer process. Assume 
we would like to simulate a fc-factor Gegenbauer process, (X t ) t , with PSD / 
as defined in (5), with Gegenbauer frequencies (Ax, . . . , A&) and long memory- 
parameters (di, . . . , dk)- The length of the realization will be N = 2 J . 

We define the band-pass variance /3? in the frequency interval Ij = [t^tt, |jtt] 
by: 

f p+l. 

0}* = 2 / 7 /(A)dA (8) 

23 + 1 

As in [7] and [2], we assume that the PSD in each frequency interval Jj, for 
which the couple (j,p) is a leaf of the tree T associated to the basis £>, is 
constant and equal to cr| . Then, the band-pass variance is (approximately) 
equal to: 

^ = 2 /^^ = 2-^L (9) 

J 2J+T 

Thus the variance of each WP coefficient is given by V[M /r J(n)] = <rj = 2 J 75| p , 
V is the variance operator. To simulate N observations of a Gegenbauer process 
(A" t ) t= i r .. ) jv with PSD /, we use the following procedure: 
1: Given an appropriate basis B and its associated tree T, calculate the 

band-pass variances /3| , (j,p) € T as in (8); 
2: For each (j,p) G T, generate 2 J_JI realizations of Wj(n), an independent 

Gaussian random variable with zero mean and variance equal to cr|„; 
3: Organize the WP coefficients Wj(n), for (j,p) E T and n = 1, . . . , 2 J ', in 
a vector Wg, and apply the the inverse DWPT to obtain the observation 
vector X= (X u . . . , X N ) T . 

In the following subsection, we examine the construction of what we term an 
appropriate basis B. 

3.2 Best-basis construction algorithm 

3.2.1 Approximate Diagonalization in a B est- Ortho- basis 

Let (Xt)t be a stationary Gegenbauer process and V its covariance matrix. Let 
7ij [£>] the entries of T [B] ; the covariance matrix of the coordinates Wg of 
(X t ) t in the ortho-basis B. One can define diagonalization as an optimization 
of the functional [17]: 

max£(£) = max5>( 7 ^D (10) 
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where e is taken as a strictly convex cost function. In practice, the optimiza- 
tion formulation of diagonalization is not widely used, presumably because it 
generally does not help in computing diagonalizations. Optimization of an ar- 
bitrary objective S over finite libraries of orthogonal bases - the cosine packets 
library and the wavelet packets library - is not a problem with good algorith- 
mic solutions. Wickerhauser [27] suggested applying these libraries in problems 
related to covariance estimation. He proposed the notion of selecting a "best 
basis" for representing a covariance by optimization of the "entropy func- 
tional" e#(7) = — log 7 over all bases in a restricted library. Authors in [16], 
developed a proposal which uses the specific choice 62(7) = 7 2 . 

In the Wickerhauser formulation, one is optimizing over a finite library and 
there will not generally be a basis in this library which exactly diagonalizes 
T. Then different strictly convex functions e(7) may end up picking different 
bases. For example, the quadratic cost function has a special interpretation 
in this context as it leads to a basis which best diagonalizes T in a least- 
squares sense [17], and is closely related to the Hilbert-Schmidt (HS) norm of 
the diagonalization error. Similarly, the -log "entropy functional" is connected 
to the Kullback-Leibler divergence [17]. Even if the approach developed in 
[16,17] was specialized to the case of 62, it is not really tied to the specific 
entropy measure; other additive convex measures can be accommodated such 
as the l a norm a > 2 or the neg-entropy, and the CW proposal makes equally 
sense. This was the starting point of our work. 



3.2.2 Proposed Algorithm 

The optimization problem of S over bases can be re-expressed as an opti- 
mization over trees, as follows. Set <^V[Wj'] = E nj e V Then as 
Euee = E(j»er En,, one is actually trying to optimize: 

E *W\ (11) 

U,p)eT 

over all recursive dyadic partitions of the spectral axis. The best basis B is 
then the one that maximizes some measure of the wavelet packets variances, 
among all the bases that can be constructed from the tree-structured library. 
The construction of the best basis can be accomplished efficiently using the 
recursive bottom-up CW algorithm defined by [1]: 

f^u^t 1 if Mw%i} + Mw%t 1 }>Mw?), 
3 m if Mw%i] + MwfS 1 } < 1 ) 



The chosen criterion lies on the comparison between some measure of WP 
coefficients variances at the children nodes and their parents. Beside the fact 
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that these variances, and then the basis, will depend on the long-memory pa- 
rameter (or even the wavelet), there is another even more important reason 
that prevents from a crude use of such a search algorithm with the cost func- 
tions that we defined above (such as e-z). Indeed, the band-pass variance of 
any node is equal to the sum of those of its children. Hence, it is not a difficult 
matter to check that any strictly convex cost functional such as those speci- 
fied above, e.g. €2 or e#, will systematically provide the basis corresponding 
to the finest partition of the spectral axis, which is clearly the worst in terms 
of complexity (i.e. number of wavelet packets). Again, this makes the wavelet 
packets machinery only of limited interest. 

Therefore, motivated by the above discussion, we were led to define, for the 
wavelet packets Wj+i and Wj+l , a new type of WP variance cost measure as 
follows: 

where A is a fixed positive constant (its value will depend for instance on the 
singularity frequency and will be given in the proof of Proposition 1). 

In the following, when we write (with a slight abuse of notation) that ^[W^+i] = 

or SvlWfl^ 1 } = 0, it will mean respectively that there exists a constant 

A < 1 such that V[WjJi] < AoVlWflt 1 } °r V^/ff 1 ] < A Y[W^ t ]. In these 
cases we will also use respectively the notations, 



y[w^] < viwgi 1 } and viwftt 1 } < v[wj 



2p 1 

+lJ 



Using the criterion defined above, algorithm (12) becomes 2 : 

BP= ,Bf +1 UBf + + 1 \ if ^[W^] = or MW%?} = 0, 
' B V j, otherwise. 

In the following, we use this algorithm to build the best-ortho-basis for a 
Gegenbauer process. 



3.3 The 1- factor case 



It is natural to build the best basis according to the shape of the PSD of our 
process. More precisely, the basis is a function of the location of the singu- 

2 Strictly speaking, this is no longer a CW algorithm. 
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larities. It means that in the case of 1-factor Gegenbauer process, the basis 
depends directly on the value of the Gegenbauer frequency. Using the nota- 
tions defined in the previous section, the recursive construction is summarized 
in the following proposition. 

Proposition 1 If (X t ) t is a stationary 1-factor Gegenbauer process, with pa- 
rameters (d,v,a) then, at node {j,p), if the frequency v is in the interval 
Ij = i§fi %r[? then: 

<*[W£i] = or Mw^f] = o, 
and consequently for algorithm (15): 

= Bf +1 U Bf + \\ 

Furthermore, if the frequency v is in the closure of the intervals and Ij+i, 
then: 

Sy[W^\ = and ^[Wflf] = 0, 
and consequently for algorithm (15): 



Proof: See Appendix A. 

To construct the best-ortho-basis of a 1-factor Gegenbauer process, we propose 
Algorithm 1 which proceeds according to Proposition 1 and the aggregation 
relation defined in (15). This algorithm is decomposed into two mains loops. 
The first one builds a family where the best-ortho-basis is included. The second 
loop is a pruning of the family to obtain the best-ortho-basis. This second loop 
corresponds to the second part of Proposition 1. 

The algorithm we propose is very fast involving only simple comparisons, and 
it does not require the calculation of variances of WP coefficients. To illustrate 
the computational speed of our algorithm we provide in Fig.l some computa- 
tion times to build bases using our method and the method of [2] 3 . In this 
example, we are only interested in the time needed to build the basis. These 
bases are built to simulate Gegenbauer process with a singularity located at 
1/12 and length equal to 2 J , with J = 6, . . . , 13. The solid line corresponds 
to the computation time of the algorithm we propose. The symbols '+', 'x' 
and '.' correspond to the computation time using the method of [2] in the case 
of respectively 'dblO' (Daubechies wavelet with q = 10 vanishing moments), 
'symlO' (Symmlet q = 10) and 'coif 5' (Coiflet q = 10). In every case the 

3 The experiments were run under the R environment on a 2.4GHz PC with 512MB 
RAM 
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Algorithm 1. 1-factor Best-Basis Search Algorithm 
Require: A Gegenbauer frequency v and sample size N = 2 J , 

Initialization 

for j — 0, . . . , J and p = 0, . . . , 2° — 1 do 

Tree(j,p) = 0. 
end for 

Main Loop 
for j — 1, . . . , J do 

for p = 0, 2, . . . , 2* - 2 do 

if v G [p/2i + \ (p+ l)/2-? +1 ] then 

Tree(j,p + 1) = 1 
end if 

if v G [(p + l)/2^' +1 , (p + 2) /2-'' +1 ] then 

Tree(j,p) = 1 
end if 
end for 
end for 

Pruning 

for j — 1, . . . , J do 

for p = 0, 2, . . . , 2 J - 2 do 

if z/ G [p/2 J+1 ,(p+ l)/2- ?+1 ] and max r=lr . )J _ i _ 1 . s=0i .,. i2 r_i Tree(j + 
r, 2 r p + s) > then 

Tree(j,p) = 
end if 
end for 
end for 

computation time increases with the length of the process. However this time 
increases always much faster for [2] than for the current method (the ratio 
of computation times is 10 to 300 times larger for the competitor method for 
series of length 64 to 8192). Typically, for a 8192-sample series, it takes 100 
ms to our algorithm to find the best basis while [2] algorithm requires 30 s. 



Examples 

We give two examples of construction of bases. Fig. 2. (a) depicts the basis built 
using the first part of Algorithm 1, to simulate a stationary Gegenbauer pro- 
cess with frequency v = 1/12. Fig.2.(b) shows the basis constructed in the case 
of a stationary Gegenbauer process with v = 0.375. The last case corresponds 
to the second situation of Proposition 1. One may remark that unlike the first 
case where the tree has at least one leaf at each scale, in this second case, 
because of the particular value of the Gegenbauer frequency, there exists scale 
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for which the tree has no leaf (see scale j = 2). For comparative purposes, 
observe that the basis provided by the approach of [2] (with a threshold 0.01) 
is highly dependent of the wavelet choice. For example in Fig.2.(c) ('db3'), 
one cannot have an idea of the singularity location. In Fig.2.(d) ('coif 5'), two 
singularties are apparent while only one is relevant. In both last cases, the 
basis is clearly over-partitioned. 

3.4 The k-factor case 

In this section we are interested in the general case: the construction of the 
appropriate basis to simulate a /c-factor Gegenbauer process. To achieve this 
goal, let us consider (Xfyt and (X^) t as respectively a (k — 1) -factor and a 
1-factor Gegenbauer processes. We denote (d±, u\, . . . , d k _i, u k _i) and (c4, Vk) 
the parameters of [X\) t and (Xf) t . Let B\ and £> 2 the best-ortho-bases of 
(X}) t and (Xf) t . We denote respectively T\ and T~2 the trees associated with 
the bases B\ and £>2- 

Let (Xt)t be a fc-factor Gegenbauer process with parameters (d\, Ui, . . . , d k , Vk). 
We denote B the appropriate basis and T the associated tree. Let B' be the 
family equal to the union of the bases B\ and B2 and let T be the associated 
tree. We are now ready to state the following, 

Proposition 2 Under the previous assumptions: 

(1) BcB' 

(2) Let (j,p) be node in the tree T such that there exists r* = 1, . . . , J — j 
and s* = 0, . . . , 2 r — 1 such that (j + r*, 2 r p + s) is also in the tree T' . 
Then: 

(j,p)&T and U + r*,2 r *p + s) e T 
Proof: See Appendix A. 

According to this last proposition, the best-ortho-basis of a /c-factor Gegen- 
bauer process may be built using k well chosen best-ortho-bases of 1-factor 
Gegenbauer processes. The steps outlined in Algorithm 2 allow to build the 
appropriate basis to simulate a /c-factor Gegenbauer process. This procedure 
lies on Algorithm 1 and results given in Proposition 2. 

Example 

Here, we give an example of construction of the best-ortho-basis for a 2- 
factor Gegenbauer process (X t ) t with Gegenbauer frequencies 1/12 and 1/24. 
Fig. 3. (a) and 3.(b) show the best-ortho-bases B\ and B2 of the processes (X]) t 
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Algorithm 2. fc-factor Best-Basis Search Algorithm 
Require: Gegenbauer frequencies and sample size N = 2 J , 

Initialization 

1: for Each Gegenbauer frequency v^i— 1, . . . , k do 

2: Construct the best-ortho-basis £>j and associated tree Tree^ using Algo- 
rithm 1. 
3: end for 

4: Tree = \J k i= {Treei (implemented using e.g. the logical OR operator under 
R or Matlab). 

Pruning 

5: for j — 1, . . . , J do 

6: for p = 0, 2, . . . , 2 j - 2 do 

7: if Tree(j,p) = 1 and max r=lr .. i j_ i _ 1 . s= o ! ..., 2 --i Tree(j + r, 2 r p + s) > 
then 

8: Tree(j,p) = 

9: end if 
10: end for 
11: end for 



and {Xf) t (see the previous section for construction of these bases). The fam- 
ily B* equal to B\ U B2 is given in Fig.3.(c). This family is not a basis, the 
intersections between its elements are not always empty, e.g. at depth j = 3, 
the elements at p = and p = 1 should not be considered as elements of 
the best-ortho-basis and must be pruned away. This is accomplished using the 
methodology developed above, and an appropriate basis for the process (X t ) t 
is obtained as represented in Fig.3.(d). 



3.5 Back to the original CW algorithm 

Our aim here is to shed light on our best-basis search algorithm by relating it to 
the original CW one. More precisely, we shall give an additive cost functional, 
which can be used within the CW algorithm, that is closely linked to our 
proposal in (15). Basically, the aggregation relation (15) can be thought of as 
a rule which at each level, enforces dyadic splitting of an interval /J if (and 
only if) the WP variance inside that interval is above a threshold, and the node 
corresponding to this interval is marked as a branch (non-terminal). Otherwise 
the interval is kept intact (marked as a leaf) and the children nodes in the 
tree are pruned away. Doing so, this procedure implicitly tries to track the 
packets that contain the singularities of the process. Hence, motivated by these 
observations, an additive variance cost functional satisfying this aggregation 
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rule can be defined as: 

MWf] = fat (fa > 5) (16) 

where 0j p is the band-pass variance as before and 5 is a strictly positive thresh- 
old. The original recursive (bottom-up) CW could then be used to minimize 
such a cost functional (termed as "Number above a threshold" functional in 
Wickerhauser book [20]). Unfortunately, the threshold remains an important 
issue to fix, and depends jointly on the singularity frequencies, the long mem- 
ory parameter, the WP level and even its location. It is therefore awkward 
to choose and control in general. To circumvent such a difficulty, a condi- 
tion involving the singularity frequencies can substitute for the thresholding 
condition in (16), that is: 

<?v[Wf] = fat (3 I = 1, . . . , k | v x e if) (17) 

From the above arguments, it turns out that minimizing the latter cost (with 
the CW algorithm) will provide us with the same basis as Algorithm 2. The 
main difference is that from a numerical standpoint, our construction algo- 
rithm is much faster and stable since there is no need to compute explicitly 
the band-pass variances, which avoids possible numerical integration problems 
(because of the PSD singularities). 



4 Analysis of decorrelation properties 



One of the approximations adopted to simulate the Gegenbauer processes us- 
ing the DWPT is that the coefficients inside each packet of the basis are 
uncorrelated. Strictly speaking, this is not true, although the expected range 
of correlation is rather weak as evidently shown by the numerical experiments 
in Fig. 4, in contrast to the long-range dependence of the process in the orig- 
inal domain. This section provides a theoretical result that establishes the 
asymptotic behavior of the covariance between WP coefficients for a 1-factor 
Gegenbauer process. 

Theorem 1 If ip has q > 1 vanishing moments with support [(N\ — N 2 + 
l)/2, (N2 — Ni + l)/2] and X{t) is a stationary 1-factor Gegenbauer process 
with Gegenbauer frequency v. Then the wavelet packet coefficients covariance 
Cov{W]l{k x ),W]l{k 2 )) decays as: 

• O (\2^ki - 2^k 2 \ 2d - 1 - R n- R ^), ifpi ^ and p 2 ^ 0, 

• O (12^ - 2^k 2 \ 2d ~ 1 ~ R ^ { ^ ) ), if pi = or p 2 = 0, 

• O (|2*fc! - 2^ 2 | M " 1 ), if Pl =p 2 = 0, 
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for all, j h j 2 , h and k 2 such that \2 jl ki - 2 j2 k 2 \ > (N* + l)(2 jl + 2 j2 ), with 
N* = max(A r 1 , N 2 ), and R p = qJ2k=oP k f or P 7^ 0, andp = (p'~ 1 p'~ 2 . . .p 1 p°) 2 
is the binary representation of p. In the last case, we note that j\ = j 2 = j. 

Proof: See Appendix B. 

This proposition generalizes the results given by [10] and [15] for the case of 
the FARIMA process. It makes an interesting statement about the order of 
correlation between well separated WP coefficients, by establishing that the 
covariance between W^{k\) and W^{k 2 ) decays exponentially over time and 
scale space. More precisely, the decay speed for pi ^ or p 2 ^ 0, depends on 
the regularity of the wavelet used, on the memory parameter of the process, 
and indirectly on the location of singularity through the frequency indices p\ 
and p 2 . However, keeping the same notations as in Proposition 1, the larger 
q, the wider the wavelet support and the fewer are the number of wavelet 
packet coefficients that satisfy the support condition \2^ki — 2^ 2 k 2 \ > (N* + 
1)(2 J1 + 2 J2 ). Thus, by choosing a wavelet with a large q, the rate of decay 
of autocovariance function increases, but over a subset of WP coefficients. 
One must then avoid inferring a stronger statement. Nonetheless, the effective 
support of a wavelet is smaller than the provided bound (see Lemma 2), and 
we expect a rapid decay in the WP coefficient's covariance for translations 
and dilations satisfying \2 jl h - 2 j2 k 2 \ > (N* + l)(2 jl + 2 j2 ). The following 
simulation study confirms these remarks. 



5 Simulation results and discussion 

5.1 Exact correlation of DWPT transformed series 

Suppose we take X as our input stationary Gegenbauer process vector, whose 
covariance matrix is V. If E is the best-ortho-basis provided by our algorithm, 
it follows that the covariance matrix of the transformed series in the WP 
domain is: 

T [B] = W^TWb (18) 

where Wg is the DWPT transform matrix operating on a vector X, whose 
columns are the basis elements of B. This equation gives the (exact) covariance 
structure for a given choice of wavelet (type, number of vanishing moments) 
and treatment of boundaries (e.g. periodic) since both are in Wg. 

Fig. 4. (a) depicts the original correlation matrix Q (resulting from T) for a 
Gegenbauer process vector (N = 64) with parameters d = 0.4 and v = 1/12. 
In Fig.4.(b)-(e) are shown the exact correlation matrices resulting from (18), 
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using respectively the Daubechies (q = 10), Symmlet (q = 10), Coiflet (q = 10) 
and Battle-Lemarie wavelets (q = 6). There is essentially no correlation within 
the packets that are far from the singularities. The most prominent correlation 
occurs within the packets close to the singularity. This effect is mainly caused 
by the support condition stated in Theorem 1 since packets near the singularity 
are those with smallest length. There is also some correlation between wavelet 
packets. A significant part of the correlation between two different packets 
seems to be concentrated along the boundaries between contiguous WP. The 
latter effect is a consequence of periodic boundary conditions. For example, the 
periodic boundary effect is higher for the Battle-Lemarie spline wavelet, whose 
support is 59 (compare to the series length of 64). But except boundary effects, 
this wavelet has a smallest between and intra-packet correlation particularly 
inside WP close to the singularity. This can be interpreted as a result of a 
sharper band-pass localization of the Battle-Lemarie filters, while the other 
wavelets have side lobes that yield more energy leaks between bands. 

To gain insight into these approximate diagonalizing capabilities of the DWPT, 
we conduct a larger scale experiment where four Gegenbauer processes with 
different frequencies and long memory parameters (d, v) are studied: three 1- 
factor with (0.4, 1/12), (0.2, 1/12), (0.3,0.016) and one 2-factor with (0.3, 1/40)- 
(0.3, 1/5). Again the influence of the wavelet on the correlation matrix result- 
ing from (18) is assessed. For comparative purposes, our bases are systemati- 
cally compared to those of [2] for each process (and wavelet as the best-basis 
of [2] depends also on the wavelet filter). We also need to consider a crite- 
rion to measure the quality of non-correlation. We here propose the Hilbert- 
Schmidt norm of the diagonalization error, which measures the sum of squares 
of the off-diagonal elements of the covariance matrix in the best-ortho-basis. 
As explained above, the method of [2] tends to over-partition the spectral axis 
yielding to too many packets. Hence, to penalize such configurations and make 
the comparison fair, we propose the following penalized criterion [17,18]: 

s{B) = \\n[B\-n Q f HS + \#{B) (19) 

where Q [B] is the correlation matrix resulting from (18) and Qq is the cor- 
relation matrix of a white noise, i.e. the identity matrix and A is a weight 
parameter balancing between the diagonalization error and the complexity of 
the tree associated to B as measured by the number of WP (leaves of 

the tree) in the basis. The value of the weight A is determined by considering 
two extreme cases. On the one hand, in the Shannon basis, we can assume that 
the decorrelation of the covariance matrix of the Gegenbauer process is perfect 
but the tree associated to this basis has too many leaves and the penalty term 
is the highest; thus S(Bs) = 2 J X. On the other hand, if one considers the 
basis £>o composed with only one leaf (i.e. the root packet Wq), there isn't any 
decorrelation of the covariance matrix. That is S(Bq) = \\Q — Qo\\%S' w hh ^ 
the correlation matrix of the Gegenbauer process whose variance-covariance 
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matrix is T. Equating the scores of these two extreme cases yields the following 
weight: 

A " 2' - 1 ■ 

Table. 1 summarizes the scores S obtained for each process as a function of 
the wavelet filter properties (type and number of vanishing moments). We 
here assumed time series of length N = 256. For details about the calculation 
of the exact autocovariance function of Gegenbauer processes and hence its 
associated covariance matrix, see [28,29,30]. These tables show that: 

• The basis provided by our algorithm is systematically better than the one 
given by [2], whatever the wavelet and process. Over-partitioning is clearly 
responsible for the bad performance of the approach in [2] . Meanwhile, the 
diagonalization error part (not shown here but will be in the next section) 
remains comparable for both bases. This means that our basis, with a re- 
duced number of packets, does not sacrify the diagonalization quality and 
yields a diagonalization error comparable to what would be obtained by 
over-partitioning. It is also worth pointing out that the approach of [2] fails 
in providing a basis for spline wavelets, and thus cannot be used in this 
case. The reason is that their best basis search algorithm strongly relies on 
a threshold on the wavelet packet filter gain, whose choice remains ad hoc. 

• For a given process, the criterion S decreases as the number of vanishing 
moments increases. This is in a very good agreement with our expectations 
as stated in Theorem 1. 

• From our experiments, we have also noticed that as the number of vanishing 
moments increases, the best basis provided by [2] tends towards the basis 
we propose. 

• For all processes, and among all tested wavelets, the Battle-Lemarie spline 
wavelet appears to provide the best score. This confirms our previous ob- 
servations. Nonetheless, the observed differences between wavelets become 
less salient at high number of vanishing moments. 

5.2 Simulation of Gegenbauer processes 

This section is devoted to the illustration of some simulation examples of 
Gegenbauer processes. The same Gegenbauer processes as in the previous 
section are considered. For each process, wavelet type and number of vanishing 
moments, M = 500 time series of length N = 256 were generated according 
to Section 3.1, using our basis and that provided by [2] method. For each 
simulated series, an unbiased estimate of the autocovariance function for the 
first N/2 lags was calculated. An average of the autocovariance function (over 
the M estimates) was then obtained and the associated correlation matrix 
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O was constructed. Finally, the HS norm of deviation between the true and 
averaged sample correlation matrices was computed: 

B(B) = \\n-Ci\\ 2 HS 

As previously, a penalized version of B(B) by the complexity of the tree as- 
sociated to B, as in (19) was also calculated (denoted B pen ) 4 . In order to 
determine which part of the score B pen is the largest contributor to the per- 
formance, and in order to not favour our best basis construction algorithm, 
both B and B pcn are displayed. The score B of the Hosking method [4], which 
is an exact simulation scheme, is also reported. The results are summarized 
in Table. 2. 

As revealed by these tables, the deviation error part B is comparable between 
the two best basis construction methods, but the penalized version differs 
significantly. This is caused by a fairly large difference in the "size" of the 
basis. Again, this backs up the statement that the method of [2] over-partitions 
the spectrum, and also agrees with the fact that in terms of performance 
our method generates as reasonable Gegenbauer processes as [2] with less 
number of packets. This also clearly provides a numerical support to our claim 
that good quality DWPT-based best-basis search, and then simulation, of 
Gegenbauer processes can be achieved without necessarily depending on the 
wavelet choice, just as it has been extensively done for 1/f processes using 
the DWT. But, one has to keep in mind that the quality of the reconstructed 
covariance structure (by assuming almost decorrelation of WP coefficients 
in the best-ortho-basis), compared to the true covariance of a Gegenbauer 
process will still depend on the wavelet. From this point of view (decorrelation 
performance), the numerical results observed for simulated data essentially 
confirm those reported in the previous subsection. 

Both the score B and its penalized version exhibit a decreasing tendency 
with increasing number of vanishing moments. This numerical evidence is a 
confirmation of the previous subsection findings and support our claims in 
Theorem 1. The Battle-Lemarie spline wavelet seems to perform the best (in 
terms of both B and -Bp en ); followed closely by the symmlets. The difference 
in performance between all wavelet types vanishes as q increases. 



4 Note that for our best-basis algorithm, and for a given process, the scores B pen is 
simply B plus a constant for all wavelets, as the penalty part in B pen only depends 
on the singularity frequencies. 



18 



6 Conclusion 



In this paper, we provided a new method to build approximate diagonal- 
izing bases for fc-factor Gegenbauer processes. Exploiting the intuitive fact 
that a wavelet packet library contains the basis where a Gegenbauer pro- 
cess could be (almost) whitened, our best-ortho-basis search algorithm was 
formulated in the case of 1-factor process and the fast search algorithm of 
Coifman-Wickerhauser was adapted to find this best basis. Using this frame- 
work, our methodology was posed in a well principled way and the uniqueness 
of the basis was guaranteed. Furthermore, unlike the approach [2], it is very 
fast (see simulations), does not depend on the wavelet choice, and is not very 
sensitive to the length of the time series. As the method construction of the 
best basis for simulation of a /c-factor Gegenbauer process relies on the 1-factor 
construction method, the same conclusions hold. 

Then, we studied the error of diagonalization in the best-ortho-basis. Towards 
this goal, we established the decay speed of the correlation between two WP 
coefficients. These results generalize the work of [10] and [15] provided in the 
case of FARIMA processes. The numerical evidence shown by our experimental 
study confirmed these theoretical findings. It has also shown that the algorithm 
introduced in the paper is appealing in that it provides good quality simulated 
Gegenbauer processes with computational simplicity and reduced complexity 
bases independently of the wavelet, which is a clear improvement over the 
existing method in [2]. Owing to these appealing theoretical and empirical 
properties, and given its practical simplicity, we feel the general practitioner 
will be attracted to our simulator. 

This new method of simulating Gegenbauer processes gives a new perspective 
for analyzing processes whose PSD singularities occur at any frequency in 
the Nyquist interval. In such a task, one could have the basis by knowing the 
process parameters (v in particular). Thus, our method has a direct application 
for bootstrap-based inference in the presence of Gegenbauer noise. 

A remaining important open problem is how we could extend this work if the 
question of interest becomes that of estimating the parameters of a /c-factor 
Gegenbauer given one or more sample paths of this process. This estimation 
problem can be accomplished in a maximum likelihood framework once the 
diagonalizing basis is found. In this case, the best-ortho-basis cannot be found 
by a naive straightforward application of Algorithm 2. Nevertheless, we have 
some promising directions that are now under investigation. Establishing the 
asymptotic behavior of such estimators also remains an open problem. One 
could also refine the estimation process by handling the residual correlation 
structure of the WP coefficients via explicit modeling by a low-order autore- 
gressive process as recently suggested in [31] for 1/f fractionally-differenced 
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processes. Additional research is still required and our current work is focusing 
on these directions. 



Appendix A 



Proof Proposition 1 : 



• Let's consider the node (j,p). We compute the variance of the WP coefficients 
at its two children: (j + 1, 2p) and (j + 1, 2p + 1). Without loss of generality, we 
assume that the frequency v is in the interval = [tJtt; ^tt[- Then a good 
approximation of the variance of the WP coefficient is given by the integral over 
the interval of the PSD. On this interval, a very good approximation to the 

PSD of the process /(A) = f^|2(cos2vrA - cos2vr^)|- M is given by C |A - v\~ 2d 
with Co a positive constant. 

Two different cases are then distinguished with associated values of A§: 
* Case v < 



YiW^} = C / 2p \v- \\~ 2d d\ + C 23+1 |A - v\- 2d d\ 

2J+ 1 U 



2 \ -d 



°o If.. 2 P\((_. 2 P\\ | / 2p+l J\ ( (2p + l 



1 - 2d \ \ 2i +l ) \ \ 2i +1 J J \ 2i +1 ) \ \ 2-J+ 1 

L _ u _ i \ 1 - 2i ) , where „ _ _ „ 



1 - 2d \ \ 2i+ l u 7 ] ' 2J+ 1 



C - 2H ( 2d 



2i+ l \ 2i+ 1 u / 
Co _._ 2d 2d \ Co ^2p+l 



2J+ 1 



where the last inequality is a consequence of the fact that u < 2~ ( - j+1 ^ in this 
case. 

To compute the variance of Wj*l_± we denote A* the location of the maxima 

of the PSD / over the interval Z^+i ■ As / is a non-increasing function over 

[ 2j y~ , anr]) it follows that this variance is bounded by a rectangle area (to a 
good approximation in this case): 

r 2 



YIW^ 1 } < 2^+i|2(co S 27rA* - cos2^)|- M . 
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Using the same approximation of the PSD as previously, we obtain: 

-2d 



r 2p +ii ^ C (2p + l V 



Thus 

Therefore, in this case we can write that VfTU^j 1 ] < V[W^J, and following the 
criterion defined in section 3.1 we have <£v[W-;£j~ ] = 0. Consequently, at the node 
(j,p), the algorithm (15) gives us: 

Bf = ^ 1 U^+ 1 . (20) 
* Case v > ^jk'- Using the same steps as in the first case, we prove that: 

When ^-pr ~ the rectangle approximation is no longer valid. But a good 
approximation of the PSD in the interval [^r, ffrr] can be Co|A — u\, where the 
constant Cq is the same as previously. Then, after some manipulations: 



1 3+1 1 1 - 2d \ 2J +l 



/ i 
1+ I 1 



V 



As by assumption -S^r ~ v, we have that < l|±I — y < — anc l then, 



1— 2d 

< 0. 



2.'H 

Finally, combining equations (21) and (22), we obtain VfW^ 1 ] ~ AiVfW^J, 
where: 

l-2d 

1 + |1 "^IfR 
i 



* = ^ " < 1. (23) 

1-11 



In this case we can write that V[Wj-^~ ] <C VfTU^-J, and using the criterion de- 
fined in section 3.1 we obtain A'fJUj 2 ^ 1 ] = 0. Finally, at the node (j,p), algorithm 
(15) gives us: 

B^ = Bf +1 UBf + f. (24) 
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In the case where the frequency v is in the closure of the intervals Ij+i and Ij+i i 

we have no relationship as ^VfW^-J = or A^Wj+i" 1 ] = 0, and one could not 
conclude. However, fortunately, at the depth j + 2, we still have: 

MW-Zt 1 } = o an d Mw-lf] = o. 

Then we easily obtain that for algorithm (15): 

KP _ K4p g 4p+l g 4p+2 „4p+3 
IJj — Dj +2 u D j+2 u D j+2 u a j+2 ■ 



Proof Proposition 2: 

(1) Let (j,p) be a node. We assume that this node is not in the tree T . It means 
that this node is not in the tree T\ neither in ?2 and in terms of threshold, we 
have 

MW-^il)) = or SvlWfZfil)] = and «^[W^ 1 (2)] = or ^[W^ 1 {2)] = °- 

As the tree 1~2 is associated to the best-ortho-basis of a 1-factor Gegenbauer 
process, the fact that the node (j,p) is node in the tree T2 means that the 
frequency is not in the interval I v - = [Jj, ^jj-]. Then in this interval, the 
function |2(cos27rA — cos 2irvk)\ ~ 2dk is bounded and has a maximum at fre- 
quency A* E Ij. Then, 

2 - 2 ?+i k-1 

V[W?L] < ^-|2(cos2vrA* - cos 2irv k )\~ 2dk / 2J+ ' TT |2(cos2vrA - cos 2nUi)\~ 2di dX 

J 57+r i=l 

2 

= ^-|2(cos2vrA* -cos2^)r 2cifc V[^ 1 (l)]. (25) 

and, using the same argument, 

2 

VIW^I 1 ] < ^|2(cos2vrA* - cos 27rz, fc )r 2 ^V[Vt^t 1 ( 1 )]- 

Finally, as ^ v [Wf^(l)] = or ^[^^(l)] = 0, we have V[WjJi] > Vft^" 1 ] 
or VfWjJj » V[W^ 2 +i +1 ], which means, 

= or MWf + t l ] = 0. 
Then, = Bf+\ U £>j+l 1 and so the node (j, p) is not in the tree T. Finally, 

B c £>'. 

(2) Here (j,p) and (j'+r*, 2 r *p + s*) (for s* = 0, . . . , 2 r * - 1) are in the tree T'. We 
denote r the minimum value of r* for which there exists as (s = 0, . . . , 2 r — 1) 
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such that the node (j + r, 2 r p + s) is in the tree T . 

Then the fact that the nodes (j, p) and (j + r, 2 r p + s) are in the tree T means 
that (j,p) is in T\ or in T2 and (j+r, 2 r p + s) is in 72 or in T\ (it is important to 
remark that both (j,p) and (j + r, 2 r p + s) cannot be in T\ or in 7^). Without 
loss of generality, we assume that (j, p) is in T2 and ( j + r, 2*p + s) is in 7i . All 
the calculations made in the following remain valid if we consider that (j,p) 
is in T\ and (j + r, 2 r p + s) is in 7-2. To simplify the notations, we assume also 
that there exists a s which is even. 

We denote Wf(l) and Wf(2), for j = 0, . . . , J and p = 0, 2^ - 1, the wavelet 
packet coefficients of respectively the processes (Xl) t and {Xf) t . From these 
sub-processes, we have that for the algorithm CW, 

• for the tree T\\ 

= S (l) U ^? r +S+1 (l) 
because <fv[^ +S (l)] = or Mwf + p r +s+ \l)} = 0, 



• for the tree 7^: 



B?(2) = B P (2) 



j 

because # V [Wf^(2)] = ¥[^(2)] and «^[W^ 1 (2)] = V\W%+\2)], 

We consider the intervals lf p r +s = [^±2, ^^±1] and /j; p r +s+1 = [^f+£±i, ^^+£±2] 
As the node (j,p) is in the tree 72, and as T2 the tree of a basis, the frequency 
v is not in the interval I^+ r +s U . 

We denote A* the location of the maximum of |2(cos27rA — cos 2nvk)\~ 2dk in 
the interval 7j +T ,_ P 1 + ' s ^ 2 (Note that the maximum is bounded). From (25), we 
have: 

V[<n < ^|2(cos2 ff A* -cos2^ fc )|-^ v[ .y a(1)]> 

■Z7T 



Then, as Sv[wf + P r +S (I)} = or ^v[^ +s+1 (l)] = 0, we have Y[wf + P+S ] » 
V[Wjl P+s+1 ] or V[Wf^ +s ] » Y[wf + P+S+1 ], which means, 

= or MWf + p r +s+l ] = 0. 
Then, because of the particular choice of r, 

B p = 2 \jB^ 

i=0 



Finally, the node (j,p) is not in the tree 7". 

However, the fact that the node (j + r, 2 r p + s) is in the tree T\ means that, 

Ar[<;;r s (i)] = <:;r s (i) a nd ^[w^^wi = <^ +2s+i u)- 
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As the frequency is not in the interval L 2 _^~ s : we obtain easily that 

and 



2 r+1 p+2si 
j+r+1 J 



w 2r+lp+2s 

vv j+r+1 



mi 2 r + 1 p+2s+U 
&V\W j+r+1 \ 



w 2 r +tp+2s+l 
VV j+r+1 



Finally, the node (j + r, 2 r p + s) is in the tree T . Using the argument, we show 
that the node (j + r*,2 r p + s*) is in the tree T . 



Appendix B 



To prove Theorem 1, the following preliminary lemmas are needed. 

Lemma 1 Let ip be a wavelet with q vanishing moments, and the associated high- 
pass QMF filters h and g. Then, for all j and p = 0, . . . , 2 3 — 1, the moments of the 
WP function tpp are such that: 



Mj, p (r) = [ ftf(t)dt = 5(r)S(p), 



for < r < R p 



where Ro = 1 and R p = q^jLJ^Pk f or P 0> an d P = {Pj-iPj-2 ■ ■ - PiPo)2 ^ s 
binary representation of p. 



Proof Lemma 1: Note that for p = 1 (wavelet basis), our result specializes to the 
traditional relation Mj t \(r) = for < r < q. The lemma can be proved either by 
induction in the original domain, or using explicit proof in the Fourier domain. We 
shall proceed according to the latter. By iterating the actions of the QMF filters h or 
g, from the root of the binary tree, to extract the appropriate range of frequencies, 
one can write that: 

ri-i 



Pi-k-i 



A-0 



(26) 



where the sequence of filters F Pk is chosen according to p = 2- 7 l Pj-i + 2- 7 2 Pj-2 + 
... + 2pi+p : 

h if Pfc = 
J if Pk = l 

and (p(0) / 0. 



F, 



Pk 



(27) 



For compactly supported wavelets with q vanishing moments, the associated high- 
pass filters g has q — 1 zeros at to = 0: 

^) = (l- e -) 9 P( e -) (28) 

where P(.) is a trigonometric polynomial bounded around u = 0. The number 
of vanishing moments of tpj(t) is equivalently given by the number of vanishing 
derivatives of t/Jj (u>) at u = 0, that is: 



M jlP (r 



i 



for r = 0, , 



(29) 



uj=0 



21 



• If p = 0, w^ioj) is just the product of low-pass niters, and ipj(t) = 4>j (t) the scaling 

function at depth j. Then, j\Aj^(r) = (f>j(0), which is non-zero with Rq = 1. If 
additional constraints are imposed on the wavelet choice (e.g. Coiflets), M.j,<o{r) 
might be zero for 1 < r < q. 

• If p ^ 0, from (26) we can write: 

n f 1 -^)'^) (so) 

k\p j _ k _ 1 =l 

where Q(.) is again bounded around u = 0. The number of vanishing moments is 
then given by the number of zeros at u = which is R p = qYlkPk- The lemma 
follows. 



Lemma 2 If the QMF h has a support in [A^A^], then the support of the WP 
function ipj(t) at each node (j,p) in the WP binary tree is always included in 
[-2? (N* + 1),2* (N* + 1)], with N* = max(|JVi|, |JV 2 |). 



Proof Lemma 2: This is proved by induction. We also use the fact that i/jq will 
be supported in the interval [iVi, N 2 ] and fa in [ Nl ~% 2+1 , jV2 ~^ 1+1 ] (see e.g. Mallat 
(1998) [21], Proposition 7.2). ■ 

Lemma 3 Let I be a collection of disjoint dyadic intervals I p whose union is the 
positive half line, and B = {^(t — 2^k) : < k < 2 J ~i ,I V - £ 2} is the associated 
orthonormal basis. Let h and g the QMFs as defined in (28). The vanishing mo- 
ments MjjOrn) of the inter- correlation function A^ 1 'f 2 (h) ofij^it) andip? 2 (t) € B 
satisfy: 

1) Pl ^ and p 2 ^ 0: M%%(m) = for 0<m<R Pl +R P2 , 

2) p 1 ^0oT P2 ^0: M%%(m)=0 for < m < R m ^( PuP2 ), 

3) Pl =p 2 = 0: M p j l$*(m)=0 for \<m<2q. 

Furthermore, the support ofA^Qi) is included in [- (N* + 1) (2^ + 2^) , (N* + 1) (2^ + T 2 )] 



Proof Lemma 3: By definition of the inter-correlation function, we have: 

/ ^imnt-h)dt (3i) 

As these WP functions belong to the orthonormal basis B, then at integer lags 
Af i f 2 {n) = 8{j l -j 2 )5{p l -p 2 )5(n). 

As far as the support is concerned, it is not a difficult matter to see, using Lemma 
2, that AfJ*(h) is supported in [- (N* + 1) (2^ + 2^ 2 ) , (N* + 1) (2^ + 2^ 2 )] . 

Let's now turn to the moments of A^'^ih). 

31,32 v ' 
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(1) Pl f and p 2 + 0: 

In this case, we know that ip^ and ip^ have respectively R pi and R P2 vanishing 
moments as defined in Lemma 1. Then, A^'J^/i) will have R Pl +R P2 vanishing 
moments since, 



h m kf?l {h)dh = - J j(v- uT^l (v)r/ 2 (u)dudv 

J2(-l) n ( nA j I v m - n ^l(v)dv j u n ^f 2 {u)du = ^ 

n=0 \ n J J J 



(32) 



for < m < R Pl + R P2 . Here, we used uniform convergence and continuity to 
invert the order of summation and integration. Note that the Fubini theorem 
allows us to invert the order of integrals. 

(2) pi 7^ or p2 7^ 0: Without loss of generality, assume that p\ ^ and p2 = 0. 
The same reasoning as above can be adopted to conclude that for < m < R Pl : 

J h m kf» 2 {h)dh = -J2 ("I)" (™) / v m - n ^l(v)dv J u n ^ 2 (u)du = 0. 



(33) 



(3) Pi = P2 = 0: 

In an orthonormal basis of wavelet packets, this situation is not possible unless 
ji = h = j- Thus, 



h m A^(h)dh = 2~ 3 J J h m <p (2~ j t) <f> (2~ j (t - h)) dtdh 
= 2 i( m + 1 ) f f u m (j)(vmv - u)dudv = 0, 



(34) 



for 1 < m < 2q, where the latter result is proved in [32]. 



Proof Theorem 1: Here we are interested in the covariance between the WP 
coefficients W^{k\) and W^fa)- We have: 



Cov 



E[X(t)X(sMl(t - 2^kM 2 2 {s - 2 h k 2 )dtd, 



cos(i/(i - a)) \t - s\ 2d - l tfl{t - 2 31 k 1 )^(s - 2^k 



After three changes of variables, u = t — 2 n k\ and v = s — 2 j2 k2, then u = t' and 
v = t' — h and finally a = 2 n (k\ — 2 n ~ n k2), we obtain: 



Cov 



W^{k 1 ),W^{k 2 



cos(u(h + a))\h + a\ 2d - l Af?*{h)dh. (36) 
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From Lemma 2 we know that that the support of A^'J^/i) is included in [— (2 n + 
2 j2 )(N* + 1), (2 jl + 2 j2 )(N* + 1)]. As h is in the support of AfJ* and by assumption 
a > (N* + l)(2 n + 2 32 ), we have h/a < 1. Hence, using the binomial series expansion 
of |l + — I" 1 and the fact that cos(z^(a + h)) ~ cos(z^a) for large a, it follows that: 



Cov 



~ |a| M 1 co 



t=i 



2d - 1 



a 



.71 



(37) 

We must then provide an upper bound on the integrals inside the braces. In the 
following we distinguish three different cases depending on the number of vanishing 
moments of A? 1 '? 2 according to Lemma 3, that is: 

(1) If pi ^ and p 2 / 0, then M P \%(m) = 0, for < m < R Pl + R P2 . We denote 
<7* = R pi + ii P2 . Then, using the fact that the q* first moments of A^'j^ are 
null, 

Ci|a| M - 1 -«*+i2 g .+i, 



Cov 



W£(*i),W£(fe) 



(38) 



with Ci a bounded constant, and: 



\R q * + i | = | cos (fa) | \a\ 2d 1 

. |M i (2d - 1 
< a 2d ' 



E 

oo 

E 



2d -I 

i 



-) 



i - h 



a 



^l{t)^f 2 {h)dtdh 



a 



i=q* + l 



C 2 \a\ 2d ~ l J2 P q * +l < CM 2d ~ l ~ q '~ l 



(39) 



i=i 



where = sup t h I ^— - 1 , C 2 and C3 are finite constants. Finally, 



Cov 



for \2Pki - 2^k 2 \ > {N* + 1)(2^ + T 2 ) and q* = R Pl + R 



vi ■ 



(40) 



(2) If pi ^ or p 2 / then Mp\%(m) = 0, for < m < R ma ^ PliP2 )- Following 
the same steps as above, we prove the second statement of Theorem 1 with 



R 



max(pi,p 2 )- 



(3) p\ = pi = 0: Mp\fp\(jn) = for 1 < m < 2q. In this particular case, we have 
necessarily j\ = j 2 = j. We must then upper-bound the covariance. From (37), 
we have: 

Cov [W°(h), W?(k 2 )] ~ C \a\ 2d ' 1 + Cilal 2 ^ 1 " 29 + R 2q+1 , (41) 
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where 



Co 



/ 



cos(f (/t+a))A '- (h)dh, C\ = cos(z^a) 



(2d - 1)! 




(2g)!(2d-l -2q)\ 



and 



-R2g+1 1 < | COS 




a 



2d-l- 



As previously, when a is large: 



Cq ~ cos(z/a) 



// 



<j)j(t)4>j{t - h)dtdh = 2 j cos(^a) |$(0)| 2 = 2 j cos(i/<g3) 



Finally, using a similar argument as in the previous cases, we find that for 
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Fig. 1. Computation time (in seconds) to build the best-ortho-basis. Solid line (our 
approach). The symbols '+', 'x' and '.' correspond to the computation time using 
the method of [2] in the case of respectively 'dbW' (Daubechies wavelet q = 10), 
'symlO' (Symmlet q = 10) and 'coif 5' (Coiflet q = 10). 
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Fig. 2. Best-orthos-basis for a Gegenbauer process, with v = 1/12 (first column) 
and v = 0.375 (second column), (a) Our basis v = 1/12, (b) our basis v = 1/5. 
Basis of [2] for v = 1/12 with (c) 'db3' and (e) 'coif 5' wavelets, and similarly for 
v = 0.375 (d)-(f). Black rectangles correspond to the leaves of the binary tree, and 
then to the partition of the spectral axis. 



d = 


0.4, v = 1/12, A = 20.7084 


d = 0.2, 


v = 1/12, A = 0.7428 


d = 0.3, 
1U.UOZD 


v = 0.016, A = 


di = d 2 
1/5, A = 


= 0.3, V! = 1/40,1*2 = 
6.0472 


y 


Whitcher Our basis 
basis 


Whitcher Our basis 
basis 


Whitcher 
basis 


Our basis 


Whitcher Our basis 
basis 




Daubechies 




Daubechies 




Daubechies 




Daubechies 


2 


2728.6 1494.5 


105.1 


52.3 


230.6 


141.4 


958.5 


639.1 


i 


1116.7 686.2 


47.3 


31.1 


188.0 


124.7 


511.4 


515.7 


6 


750.4 441.8 


31.9 


23.3 


141.8 


116.7 


444.1 


422.5 


8 


632.7 352.4 


28.9 


20.2 


144.0 


118.8 


408.0 


361.9 


10 


421.7 308.2 


21.9 


18.4 


140.3 


115.0 


400.5 


315.3 




Symmlet 




Symmlet 




Symmlet 




Symmlet 


4 


2211.0 677.0 


86.4 


30.8 


214.4 


120.9 


505.3 


513.8 


6 


1371.0 444.7 


54.2 


23.4 


210.4 


116.2 


445.0 


423.5 


8 


980.4 341.7 


39.1 


20.1 


178.6 


114.7 


409.0 


359.7 


10 


693.1 297.3 


28.2 


18.2 


138.8 


113.5 


401.2 


313.7 




Coiflet 




Coiflet 




Coiflet 




Coiflet 


2 


2697.3 1081.1 


106.0 


42.3 


224.8 


132.7 


713.5 


589.5 


4 


1449.4 638.4 


58.9 


29.6 


215.2 


121.0 


467.1 


500.7 


6 


841.1 412.7 


35.0 


22.4 


141.2 


116.1 


416.1 


406.2 


8 


703.2 327.8 


29.0 


19.5 


140.1 


114.9 


363.3 


342.7 


10 


581.2 287.7 


26.6 


17.6 


138.6 


113.2 


349.9 


297.8 




Battle-Lemarie 




Battle-Lemarie 


Battle-Lemarie 




Battle-Lemarie 


2 


657.3 




29.4 




121.3 




374.6 


4 


267.5 




16.1 




111.8 




238.1 


6 


247.8 




14.4 




110.7 




191.1 



Table 1 



S score as a function of the number of vanishing moments for each wavelet family. 
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Fig. 3. (a) Best-ortho-basis B\. (b) Best-ortho-basis B>2- (c) Union of bases B\ and 
£?2- (d) Best-ortho-basis for the two factor Gegenbauer process {X t )t- 

P3 (b) 




Fig. 4. (a) Correlation matrix of a Gegenbauer process with parameters d = 0.4 
and v = 1/12. Correlation matrix f2 [B] of the WP coefficients in our best-ortho-basis 
for 'dbW filter (b), 'symW filter (c), 'coif 5' filter (d), and 'bat& filter (e). 
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Table 2 

Squared difference (B) and its penalized version (B pen ) for our basis (superscript 
CF) and Whitcher basis [2] (superscript W) as a function of the number of vanishing 
moments for each wavelet family. 



